library(AIUQ)

Example 1: Simulate BM and get estimated parameters using BM model

set.seed(1)

## Simulation
sim_bm = simulation()
show(sim_bm)
## Frame size:  200 200 
## Number of time steps:  200 
## Number of particles:  50 
## Stochastic process:  BM 
## Variance of background noise:  20 
## sigma_bm:  1
## Plot simulated particle trajectory
plot_traj(sim_bm)

## Plot intensity profile for different frames 
plot_intensity(sim_bm@intensity, sz=sim_bm@sz) #first frame 
plot_intensity(sim_bm@intensity, sz=sim_bm@sz,frame=10, color=T) #10th frame, color image

## AIUQ method: use BM as fitted model 
sam = SAM(sim_object=sim_bm, model_name='BM')
show(sam)
## Fitted model:  BM 
## Number of q ring:  99 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 
## True parameters in the model:  2 
## Estimated parameters in the model:  2.010371 
## True variance of background noise:  20 
## Estimated variance of background noise:  19.90164 
## Maximum log likelihood value:  -18150726 
## Akaike information criterion score:  36301654
## Plot shifted intensity in reciprocal space 
plot_I_q = matrix(sam@I_q[,1], sam@sz[1],sam@sz[2]) #first frame 
plot3D::image2D(abs(fftshift(plot_I_q)),main="Foruier tansformed intensity")

## User can select q range via AIUQ_thr or index_q
sam = SAM(sim_object=sim_bm, AIUQ_thr=c(0.99,0.6)) 
#Note: Default model_name is "BM", it's ok to not specify this argument if want to fit with BM model 
show(sam)
## Fitted model:  BM 
## Number of q ring:  99 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 
## True parameters in the model:  2 
## Estimated parameters in the model:  2.008843 
## True variance of background noise:  20 
## Estimated variance of background noise:  20.05226 
## Maximum log likelihood value:  -8065526 
## Akaike information criterion score:  16131176
sam = SAM(sim_object=sim_bm, index_q_AIUQ=5:50)
show(sam)
## Fitted model:  BM 
## Number of q ring:  99 
## Index of wave number selected:  5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 
## True parameters in the model:  2 
## Estimated parameters in the model:  2.015093 
## True variance of background noise:  20 
## Estimated variance of background noise:  20.02955 
## Maximum log likelihood value:  -6198715 
## Akaike information criterion score:  12397526

Example 2: Simulate OU process and get estimated parameters using OU model

set.seed(1)

## Simulation
sim_ou = simulation(sigma_ou=4, model_name="OU")
show(sim_ou)
## Frame size:  200 200 
## Number of time steps:  200 
## Number of particles:  50 
## Stochastic process:  OU 
## Variance of background noise:  20 
## (rho, sigma_ou):  0.95 4
## Plot simulated particle trajectory
plot_traj(sim_ou)

## AIUQ method: fitting using OU model
sam_ou = SAM(sim_object=sim_ou, model_name=sim_ou@model_name) 
show(sam_ou)
## Fitted model:  OU 
## Number of q ring:  99 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 
## True parameters in the model:  0.95 64 
## Estimated parameters in the model:  0.9499009 64.42877 
## True variance of background noise:  20 
## Estimated variance of background noise:  19.84521 
## Maximum log likelihood value:  -18293754 
## Akaike information criterion score:  36587711
## Plot true MSD and estimated MSD with uncertainty
plot_MSD(object=sam_ou, msd_truth=sam_ou@msd_truth) #in log10 scale
plot_MSD(object=sam_ou, msd_truth=sam_ou@msd_truth,log10=F) #in real scale

Example 3: Simulate BM with smaller frame size and length and get estimated parameters with uncertainty using BM model

set.seed(1)

## Simulation
sim_bm = simulation(sz=100,len_t=100,sigma_bm=0.5)
show(sim_bm)
## Frame size:  100 100 
## Number of time steps:  100 
## Number of particles:  50 
## Stochastic process:  BM 
## Variance of background noise:  20 
## sigma_bm:  0.5
## Plot simulated particle trajectory
plot_traj(sim_bm)

## AIUQ method: fitting using BM model
sam = SAM(sim_object=sim_bm, uncertainty=T)
show(sam)
## Fitted model:  BM 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  0.5 
## Estimated parameters in the model:  0.5050874 
## True variance of background noise:  20 
## Estimated variance of background noise:  20.08774 
## Maximum log likelihood value:  -2289361 
## Akaike information criterion score:  4578825
## Plot true MSD and estimated MSD with uncertainty 
plot_MSD(object=sam, msd_truth=sam@msd_truth) #in log10 scale 
plot_MSD(object=sam, msd_truth=sam@msd_truth,log10=F) #in real scale

Example 4: Load csv file with intensity profile being a matrix with dim time by space*space (T_SS_mat) and est parameters using BM model

## Load intensity profile from csv file 
file_loc = "/Users/yue/Downloads/"
file_name = "intensity_record_BM_beta_0_02_B_40_len_100_frame_size_100"
intensity = data.table::fread(paste(file_loc,file_name,".csv",sep=""))
intensity = intensity[,-1]

## AIUQ method: use BM as fitted model 
mindt = 1 #minimum lag time
pxsz = 1 #pixel size
sz=c(100,100) #frame size of the image 
sam = SAM(intensity=intensity,mindt=mindt,pxsz=pxsz,sz=sz)
## start of another optimization, initial values:  0.3625398 3.5165
show(sam)
## Fitted model:  BM 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  NA 
## Estimated parameters in the model:  0.02040957 
## True variance of background noise:  NA 
## Estimated variance of background noise:  20.11201 
## Maximum log likelihood value:  -2102329 
## Akaike information criterion score:  4204761
#sam@msd_est #estimated MSD

## DDM method: use fixed A and B  
sam = SAM(intensity=intensity,mindt=mindt,pxsz=pxsz,sz=sz, method="DDM_fixedAB")
show(sam)
## Fitted model:  BM 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  NA 
## Estimated parameters in the model:  8.336566 
## True variance of background noise:  NA 
## Estimated variance of background noise:  20.35795

Example 5: Load tif file (SST_array) from real experiment and est parameters using FBM model

## Load real experiment video 
file_loc = "/Users/yue/Desktop/Research/2.SAM/AIUQ_package/data/"
file_name = "4pct_PVA_100nm_25C"
intensity = bioimagetools::readTIF(paste(file_loc,file_name,".tif",sep=""),as.is=TRUE)

## Plot intensity profile for different frames 
intensity_str = 'SST_array' #intensity profile structure
plot_intensity(intensity, intensity_str=intensity_str) #first frame 
plot_intensity(intensity, intensity_str=intensity_str, frame=15) #15th frame

## AIUQ method: fitting using FBM model 
mindt = 0.0309 #minimum lag time
pxsz = 0.2930  #pixel size
model_name = "FBM"
sam = SAM(intensity=intensity,intensity_str=intensity_str,mindt=mindt, pxsz=pxsz, uncertainty=T, AIUQ_thr=c(0.99,0), model_name=model_name)
show(sam)
## Fitted model:  FBM 
## Number of q ring:  255 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 
## True parameters in the model:  NA 
## Estimated parameters in the model:  0.634021 1.014347 
## True variance of background noise:  NA 
## Estimated variance of background noise:  150.4655 
## Maximum log likelihood value:  -41045311 
## Akaike information criterion score:  82090786
## Plot reference MSD vs estimated MSD with uncertainty 
msd_truth = 0.7589*sam@d_input
plot_MSD(object=sam, msd_truth=msd_truth)

Example 6: Simulate FBM and compare estimated parameters using SAM and DDM methods

set.seed(1)
## Simulation
sim_fbm = simulation(sz=100,len_t=100,model_name="FBM")
show(sim_fbm)
## Frame size:  100 100 
## Number of time steps:  100 
## Number of particles:  50 
## Stochastic process:  FBM 
## Variance of background noise:  20 
## (sigma_fbm, Hurst parameter):  2 0.3
## AIUQ method
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name)
show(sam)
## Fitted model:  FBM 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  8 0.6 
## Estimated parameters in the model:  8.02211 0.5969453 
## True variance of background noise:  20 
## Estimated variance of background noise:  15.05404 
## Maximum log likelihood value:  -2414150 
## Akaike information criterion score:  4828404
# Plot true MSD vs estimated MSD 
plot_MSD(object=sam,msd_truth=sam@msd_truth)

## DDM method with fixed A and B, using all q 
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_fixedAB")
show(sam)
## Fitted model:  FBM 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  8 0.6 
## Estimated parameters in the model:  6.561253 0.8439145 
## True variance of background noise:  20 
## Estimated variance of background noise:  20.55181
## DDM method with fixed A and B, using user defined q range 
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_fixedAB", index_q_DDM=3:40)
show(sam)
## Fitted model:  FBM 
## Number of q ring:  49 
## Index of wave number selected:  3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 
## True parameters in the model:  8 0.6 
## Estimated parameters in the model:  7.04659 0.8025151 
## True variance of background noise:  20 
## Estimated variance of background noise:  20.55181
## Add line of estimated MSD from DDM method with fixed A and B and user defined q range
lines(log10(sam@d_input[-1]),log10(sam@msd_est[-1]),type='p',
              col='red', pch=24, cex=0.8)

## DDM method with estimated A and B,using all q
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_estAB")
show(sam)
## Fitted model:  FBM 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  8 0.6 
## Estimated parameters in the model:  4.001553 0.9926345 
## True variance of background noise:  20 
## Estimated variance of background noise:  37.45107
## DDM method with estimated A and, using user defined q range 
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_estAB", index_q_DDM=3:40)
show(sam)
## Fitted model:  FBM 
## Number of q ring:  49 
## Index of wave number selected:  3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 
## True parameters in the model:  8 0.6 
## Estimated parameters in the model:  5.036552 0.9891863 
## True variance of background noise:  20 
## Estimated variance of background noise:  37.86241

Example 7: User defined MSD structure

set.seed(1)

## Simulation
sim_bm = simulation(sz=100,len_t=100,sigma_bm=0.5)
show(sim_bm)
## Frame size:  100 100 
## Number of time steps:  100 
## Number of particles:  50 
## Stochastic process:  BM 
## Variance of background noise:  20 
## sigma_bm:  0.5
## User defined MSD structure: function of parameters and 
#                                          vector of lag times 
msd_fn = function(param, d_input){
  MSD = param[1]*d_input+param[2]*d_input^2
}

# show MSD and MSD gradient with a simple example 
theta = c(2,1)
d_input = 0:10
model_name = "user_defined"
MSD_list = get_MSD_with_grad(theta=theta,d_input=d_input,model_name=model_name,
                             msd_fn=msd_fn)
MSD_list$msd
##  [1]   0   3   8  15  24  35  48  63  80  99 120
## AIUQ method: fitting using user_defined model 
sam = SAM(sim_object=sim_bm, model_name=model_name, msd_fn=msd_fn, num_param=2)
show(sam)
## Fitted model:  user_defined 
## Number of q ring:  49 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
## True parameters in the model:  0.5 
## Estimated parameters in the model:  0.5047796 0.000272548 
## True variance of background noise:  20 
## Estimated variance of background noise:  20.08427 
## Maximum log likelihood value:  -2289361 
## Akaike information criterion score:  4578827

Example 8: Extract empirical/modeled Dqt and ISF outside/within SAM

set.seed(1)

## Simulation
sim_bm = simulation()
show(sim_bm)
## Frame size:  200 200 
## Number of time steps:  200 
## Number of particles:  50 
## Stochastic process:  BM 
## Variance of background noise:  20 
## sigma_bm:  1
## AIUQ method: use BM as fitted model 
sam = SAM(sim_object=sim_bm, model_name='BM')
show(sam)
## Fitted model:  BM 
## Number of q ring:  99 
## Index of wave number selected:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 
## True parameters in the model:  2 
## Estimated parameters in the model:  2.010371 
## True variance of background noise:  20 
## Estimated variance of background noise:  19.90164 
## Maximum log likelihood value:  -18150726 
## Akaike information criterion score:  36301654
# Compute empirical/modeled Dqt and ISF after SAM class is constructed 
par(mfrow = c(2, 2))
## Compute observed dynamic image structure function (Dqt)
Dqt = get_dqt(sam) #q index range is all qs without user specify 
Dqt = get_dqt(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6, main="Observed dynamic image structure function")

## Compute empirical intermediate scattering function (ISF)
ISF = get_isf(sam) #q index range is all qs without user specify 
ISF = get_isf(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Empirical intermediate scattering function")

## Compute modeled Dqt 
modeled_Dqt = modeled_dqt(sam) #q index range is all qs without user specify 
modeled_Dqt = modeled_dqt(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(modeled_Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6, main="Modeled dynamic image structure function")

## Compute modeled ISF 
modeled_ISF = modeled_isf(sam) #q index range is all qs without user specify 
modeled_ISF = modeled_isf(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(modeled_ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Modeled intermediate scattering function")

# Compute empirical/modeled within SAM class

## AIUQ method: use BM as fitted model and compute empirical Dqt within class  
sam = SAM(sim_object=sim_bm, model_name='BM', output_dqt=T)
matplot(sam@d_input[-1], t(sam@Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6,main="Observed dynamic image structure function")

## AIUQ method: use BM as fitted model and compute empirical Dqt,ISF within class 
sam = SAM(sim_object=sim_bm, model_name='BM', output_isf=T) # Note: empirical ISF is calculated through observed Dqt. By setting output_isf is true, Dqt is also computed even without "output_dqt=T".
matplot(sam@d_input[-1], t(sam@Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6,main="Observed dynamic image structure function")
matplot(sam@d_input[-1], t(sam@ISF[-c(51:99),]), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Empirical intermediate scattering function")
matplot(sam@d_input[-1], t(sam@ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Empirical intermediate scattering function")

## AIUQ method: use BM as fitted model and compute ISF within class  
sam = SAM(sim_object=sim_bm, model_name='BM', output_modeled_isf=T)
matplot(sam@d_input[-1], t(sam@modeled_ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6,main="Modeled intermediate scattering function")

## AIUQ method: use BM as fitted model and modeled Dqt,ISF within class 
sam = SAM(sim_object=sim_bm, model_name='BM', output_modeled_dqt=T) # Note: Modeled Dqt is calculated through modeled ISF. By setting output_modeled_dqt is true, modeled_ISF is also computed even without "output_modeled_isf=T".
matplot(sam@d_input[-1], t(sam@modeled_ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Modeled intermediate scattering function")
matplot(sam@d_input[-1], t(sam@modeled_Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6,main="Modeled dynamic image structure function")